type vector = {mutable x: float;mutable y: float;mutable z: float}
type point = vector
type triangle = {sa: point; sb: point; sc: point}
type trianglei = {isa: int; isb: int; isc: int}


(* cross
 * 3D cross product : u x v = w
 *)
let cross (u: vector) (v: vector) =
	{x = u.y*.v.z -. u.z*.v.y;
	y = u.z*.v.x -. u.x*.v.z;
	z = u.x*.v.y -. v.x*.u.y}

let magnitude (u: vector) =
  sqrt ((u.x*.u.x) +.
        (u.y*.u.y) +.
        (u.z*.u.z))

(* normalize
 * Normalize a vector : |u| = 1
 *)
let normalize (u: vector) =
	let magn = sqrt ((u.x*.u.x) +.
			 (u.y*.u.y) +.
			 (u.z*.u.z))
	in
	u.x <- u.x /. magn;
	u.y <- u.y /. magn;
	u.z <- u.z /. magn

let normalized (u: vector) =
	let magn = sqrt ((u.x*.u.x) +.
			 (u.y*.u.y) +.
			 (u.z*.u.z))
	in
          { x = u.x /. magn;
	y = u.y /. magn;
	z = u.z /. magn
          }
(* fromC3
* Converts a triple float into a vector
*)
let fromC3 (vx,vy,vz) =
	{x = vx; y = vy; z = vz}
(* toC3
* Converts a vector into a triple float
*)
let toC3 (v: vector) =
	(v.x, v.y, v.z)

(* diff
* computes that vector substraction u - v
*)

let diff (u: vector) (v: vector) =
	{x = u.x -. v.x; y = u.y -. v.y; z = u.z -. v.z}

let add (u: vector) (v: vector) =
	{x = u.x +. v.x; y = u.y +. v.y; z = u.z +. v.z}
let mult n (u: vector) =
	{x = u.x *.n; y = u.y *.n; z = u.z *.n}
(* dot2d
 * 2d dot product:  vector(ua,ub) * vector(va,vb)
 *)
let dot2d (ua: point) (ub: point) (va: point) (vb: point) =
        (ub.x-.ua.x)*.(vb.x-.va.x) +. (ub.y-.ua.y)*.(vb.y-.va.y)


(* is_left
 * Determines the position of a point according to an
 * infinite line using perpendicular product.
 * Inputs:
         * a as point : first point in the line
         * b as point : second point in the line
         * c as point : point to be tested
 * Output:
         * float : val > 0  => point is in the left halfspace
                   val < 0  => point is in the right halfspace
                   val = 0  => point is in the line
 *)
let is_left (a: point) (b: point) (c: point) =
        (b.x -. a.x)*.(c.y -. a.y) -. (c.x -. a.x)*.(b.y-.a.y)


(* pt_in_seg
 * Tests if one point is  contained by a segment,
 * using square dot product.
 * Inputs:
         * a as point : first point of the segment
         * b as point : second point of the segment
         * c as point : point to be tested
 * Output:
         * boolean : true  => point is in segment
                     false => point is out of the segment
 *)
let  pt_in_seg (a: point) (b: point) (c: point) =
        let d = dot2d a b a c in
        let abab = dot2d a b a b in
           if  d>=0. && d<= abab then
              d*.d = abab *.(dot2d a c a c)
           else
              false


(* pt_in_tborders
 * Determines if a given point is in at least on one
 * of a triangle border's.
 * Inputs:
         * t as triangle : triangle to be tested
         * p as point : point to be tested
 * Outputs:
	 * val = 0 iff pt is not in any border
	 * val = 1 iff pt is in [sasb]
	 * val = 2 iff pt is in [sbsc]
	 * val = 3 iff pt is in [scsa]
 *)
let  pt_in_tborders (t: triangle) (p: point) =
   if  pt_in_seg t.sa t.sb p then
      1
   else
      if pt_in_seg t.sb t.sc p then
     	 2
      else
         if pt_in_seg t.sc t.sa p then
   	   3
	 else
	   0

(* pti_in_tiborders
 * Determines if a given point is in at least on one
 * of a triangle border's.
 * Inputs:
         * f as int -> pointi : function transforming an
                               int index to a point.
         * t as trianglei : index-triangle to be tested
         * p as int : index of the point to be tested
 * Outputs:
	 * val = 0 iff pt is not in any border
	 * val = 1 iff pt is in [isaisb]
	 * val = 2 iff pt is in [isbisc]
	 * val = 3 iff pt is in [iscisa]
 *)
let pti_in_tiborders f t p =
   pt_in_tborders {sa= f t.isa;sb= f t.isb; sc= f t.isc} p

(* triangle_split_seg
 * Splits a CW triangle (defined by pointe a, b ,c)
 * into two CW sub-triangles. The original
 * triangle is splited by inserting d on the segment [ab].
 * Return reliable resuls iff d is in [ab].
 * Inputs:
	 * a as point : first point of the splited triangle.
	 	        Must be the first point of the future
			splitted edge.
	 * b as point : second point of the splited triangle.
	 	        Must be the first point of the future
			splitted edge.
	 * c as point : third point of the splited triangle,
		        common to the two new ceated triangles.
	 * d as point : splitting point.
 * Outputs:
	 * val as (triangle,triangle) :
		couple of CW sub-triangles: adc, dbc
 *)
let triangle_split_seg (a: point) (b: point) (c: point) (d: point) =
	({sa=a;sb=d;sc=c} , {sa=d;sb=b;sc=c})

(* trianglei_split_segi
 * Same as triangle_split_seg with triangli indexed-vertices triangles.
 *)
let trianglei_split_segi a b c d =
	({isa=a;isb=d;isc=c} , {isa=d;isb=b;isc=c})

(* trianglei_split_segi_n
 * Splits the n-th edge of a given trianglei.
 * Uses trianglei_split_segi function to split.
 *)
let trianglei_split_segi_n t d = function
    | 1 -> trianglei_split_segi t.isa t.isb t.isc d
    | 2 -> trianglei_split_segi t.isb t.isc t.isa d
    | _ -> (* 3 *) trianglei_split_segi t.isc t.isa t.isb d

(* triangle_split
 * Splits a CW triangle to form two CW sub-triangles. The original
 * triangle is splited by inserting p.
 * Returns reliable resuls iff c is in one of the triangle's edges.
 * Inputs:
	 * t as triangle : triangle to split.
	 * p as point : splitting point.
	 * _ as int : edge identifier. Can be:
		1 => splits edge ab
		2 => splits edge bc
		3 => splits edge ca
 * Outputs:
	 * val as (triangle,triangle) : couple of CW sub-triangles.
 *)
let triangle_split (t: triangle) (p: point) = function
 	  1 -> triangle_split_seg t.sa t.sb t.sc p
	| 2 -> triangle_split_seg t.sb t.sc t.sa p
	| _ -> triangle_split_seg t.sc t.sa t.sb p

(* trianglei_split
 * Same as triangle_split with trianglei indexed-vertices triangles.
 *)
let trianglei_split t p = function
 	  1 -> trianglei_split_segi t.isa t.isb t.isc p
	| 2 -> trianglei_split_segi t.isb t.isc t.isa p
	| _ -> trianglei_split_segi t.isc t.isa t.isb p

